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Abstract 

The entanglement entropy of a subsystem A of a quantum system is expressed, in 
the replica approach, through analytic continuation with respect to n of the trace of the 
n— th power of the reduced density matrix. This trace can be thought of as the vacuum 
expectation value of a suitable observable in a system made with n independent copies 
of the original system. We use this property to numerically evaluate it in some two- 
dimensional critical systems, where it can be compared with the results of Calabrese 
and Cardy, who wrote the same quantity in terms of correlation functions of twist 
fields of a conformal field theory. Although the two calculations match perfectly even 
in finite systems when the system A consists of a single interval, they disagree whenever 
the subsystem A is composed of more than one connected part. The reasons of this 
disagreement are explained. 



1 Introduction 



In a quantum system, performing a local measurement may instantaneously affect local 
measurements far away. This is a manifestation of the quantum entanglement. A useful 
measure of this property in extended quantum systems with many degrees of freedom is the 
von Neumann or entanglement entropy: one considers a pure quantum state (typically 
the ground state) that can be subdivided into two subsystems A and B and constructs the 
reduced density matrix 

p A = ti B \m)(m\ (i.i) 

by tracing over the degrees of freedom of B. The von Neumann or entanglement entropy is 
defined to be 

S A = -tr p A lnp A = -trp B lnp B = S B • (1.2) 

Of particular interest is the case in which the two subsystems A and B correspond to 
two connected regions (the "outside" and the "inside") in the space-time, where it was 
argued that the entanglement (or geometric) entropy is deeply related to the physics of the 
black holesp], [2} [3]. An important element of the present understanding is its holographic 
interpretation [5]. This point of view seems to indicate a purely geometric way of computing 
the entanglement entropy in strongly coupled conformal field theories [6]. 

The entanglement entropy has been also extensively studied in low dimensional quan- 
tum systems as a new tool to investigate the nature of quantum criticality [H El EE UTil fTT] . 
Several different calculations based on the conformal field theory (CFT) describing the uni- 
versal properties of the quantum phase transition describing 1+1 dimensional systems, like 
quantum spin chains, have shown that the entropy grows logarithmically with the size I of 
the subsystem A as [Zl El El HD] 

S A = ~M£/a), (1.3) 

where c is the conformal anomaly and a an ultraviolet cutoff. 

The trace of the n-th power of the reduced density matrix, which could be identified with 
the Tsallis entropy [T2], plays a major role in the replica approach to entanglement entropy 
[TO] , yielding 

S A = -\im^-trp n A . (1.4) 
an 

Our goal in this paper is to describe a simple method to directly measure tr p n A in whatever 
local lattice field theory in any space-time dimension, based on the observation that this 
quantity can be expressed, as we shall explain below, as the vacuum expectation value of 
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a suitable observable defined in a system made with n independent copies of the original 
quantum system. 

We apply our method to check various consequences of the replica approach in some 
(1 + 1) -dimensional critical quantum system described by a relativistic conformal field the- 
ory (CFT), including a system with c < 0, where the entanglement entropy appears to be 
negative, finding also in that case complete agreement: evidently in those non-unitary sys- 
tems the mixed state obtained by tracing over the degrees of freedom of B is more ordered 
than the pure state |\&). 

When the subsystem A consists of N disjoint intervals A = (ui, Ui) . . . (un, Vn) the quan- 
tity tr p n A is proportional, at criticality, to the n— th power of the correlation function of 2N 
local primary operators of complex scaling dimensions 

sitting on the end points of the N intervals: 

N 

tr^oc(n^(%,%)^(^^)) n - (1.6) 

Notice that these primary operators do not belong to the Kac table. They can be considered 
as a special kind of twist fields, called branch-point twist fields [13] , because they are naturally 
related to the branch points in the n-sheeted Riemann surface where the system is defined in 
its Euclidean functional integral formulation. These primary operators also carry a conserved 
charge related to the orientation of the intervals. 

The proof of the above statements relies (among other things) on the clever observation 
[TO] that the vacuum expectation value of the stress tensor T(z) near a branch point has 
the same functional dependence as the one generated by a primary operator with scaling 
dimensions A n . Since the same branch point is present in each sheet of the Riemann surface, 
it is easy to conclude that tr p\ has the form ( 11. 6ft . 

By making a further assumption we shall discuss in the next Section, the authors of [10] 
proposed an explicit functional form0 



trp^ cx 



Uj,k<N( V k ~ Uj) 



4nA n 

(1-7) 



1 There is a typo in this formula in [TU]. We thank P.Calabrese for pointing this out to us. 
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Note that the points of type u and of type v, associated respectively with <3>~ and are 
treated in a different way. This is justified only when the intervals are intrinsically oriented, 
which is not the case when n = 2, where all the branch points are on the same footing 
and $2 ( w ) = ^t( w )- We shall discuss below the N = n = 2 case in detail. When n 
is an odd number and c > we shall argue, and corroborate with numerical experiments, 
that whenever two twist fields of the same charge approach each other the above quantity 
develops an ultraviolet divergence instead of vanishing as (11.71) would require. 

The origin of such a discrepancy can be traced to the assumption that the only singular- 
ities of T(z) are the double poles at the branch points of the multi-sheeted Riemann surface, 
which is not the case of the conformal mappings w = f(z) considered in [10J. This point will 
be further discussed in the next Section. 

The contents of this paper are as follows. In the next Section we summarise the main 
results of the replica approach to entanglement entropy in CFT and discuss some difficulties 
in the use of the conformal mappings when the number of branch points is larger than 
two. In the following Section we describe our method of evaluating trp^ as the vacuum 
expectation value of a suitable observable; as a byproduct we uncover a special quantum 
symmetry, which is exact only at zero temperature, yielding at once the important identity 
trp^ = trp^ and other useful relationships. In Section H] we implement the method in some 
specific two-dimensional systems, namely the Potts models: they can be easily simulated at 
their critical point and the corresponding CFT is known. We finish with some conclusions. 

2 Conformal maps and twist fields 

Translational and rotational invariance of a CFT on the complex z— plane C imply (T(z))c — 
0. Under a conformal mapping w — > z = f(w) the stress tensor T transforms as 



is called Schwarzian derivative. Comparing the vacuum expectation value of both sides of 





where 




(2.2) 



(T2~Tj) yields 



(2.3) 



where 71 is the Riemann surface associated with f(w). 
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The identity {/, w} = —^VT'-^^hr implies at once that {/, w} = if and only if f(w) is 
a linear fractional transformation, which is the only conformal mapping of the whole C onto 
itself. In any other case (T(w)) ^ and some symmetry of the original system is lost. In 
particular, using again the above-mentioned identity, it is easy to see that the only conformal 
mappings conserving translational invariance ( hence (T(w))-ji = const.) are the exponential 
mappings, which can be written in a suitable basis in the form 

w 

z = exp(27r-) , (2.4) 

which maps the cylinder composed of the infinite strip < Qmw = y < L with periodic 
boundary conditions into the whole z plane excepted the origin. As a parenthetic remark, we 
note that this exponential mapping is the key ingredient to find universal finite size effects 
of two-dimensional field theories at criticality [13]. In any other conformal mapping, being 
(T(w)) a non-constant analytic function, it has to be singular somewhere. By dimensional 
reasons, the isolated singularities of (T(w)) are double poles. 

As a simple illustrative example, apply (12.31) to the mapping z = f(w) = It yields 
at once 

1/n 2 ) 
24 w 2 

where TZ n is now a n— sheeted Riemann surface branched over and oo. We can now first 
displace these points with a linear fractional transformation g(w) = an d then consider 
the composed transformation 

(in — 1l\ n 
. (2.6) 
w — v J 

Being g a linear fractional transformation, we may take advantage of the remarkable identity 
{f(g),w} = g' 2 {f,g} to get 

{T{W)) ^ = ° 24 (w-unw-v) 2 ■ (2 - 7) 



As first pointed out in [10], this expression coincides, up to a normalising constant, with 
the three-point function (T(w)$> n (u, u)$> n (v, v))c where Q n (w,w) is a primary operator, 
called branch-point twist field [13], with complex scaling dimensions A n = A n = c 1 "^" 
associated with the branch points w = u and w = v of TZ n , as already anticipated in the 
Introduction. Combining such an observation with the fact that the surface constructed with 
the replica method has the same intrinsic geometric properties as 1Z n Calabrese and Cardy 
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[TO] were able to argue that 



u — V 



\(n-l/n) 



(21 



where a is an UV cutoff and A is the interval joining the branch points u and v. 

To make contact with numerical simulations it is useful to rewrite this correlator in a finite 
geometry using the transformation law of the primary operators under a suitable conformal 
mapping w = h(() : 

(*;(Ci,CO*£(Ca,C0>w= \™\(i)w\( 2 )\ 2An (K(u,u)$+(v,v))c , (2.9) 

where £1 and (2 are the inverse images of u and v and 7?. is the Riemann surface associated 
with h. Using the exponential mapping (j2.4p and arranging for instance the branch points 
to lie on the imaginary axis of the cylinder ( = x + % y, one finally finds [TO] 



sin(y/L) J ' (2 - 10) 

with £ = 1 2/2 - 2/1 1 - 

We checked this formula with numerical experiments in two different critical systems, 
one, the Ising model, with positive conformal anomaly (c = ~) and the other, a Q-state 
Potts model with Q < 1, corresponding to a conformal theory with a negative conformal 
anomaly (c = — In either case we studied the dependence on n, L and £, finding perfect 
agreement with (12.101) provided the distances of the branch points and the sizes of the lattice 
are large enough with respect to the lattice spacing (see Section H] for more details). 

2.1 General case 

When A consists of more than one disjoint interval, i.e. A = (ui, t>i) . . . (u^, v^), it is natural 
to assume that the Tsallis entropy tr p n A be proportional to the n— th power of the correlation 
function of the twist fields associated with the 2N branch points (see Eq. fll.6l) ). 

The Ward identities of the conformal invariance tell us that the only singularities of 
(T(w)$~(mi, Mi) . . . $+(t>Ar, vn)) are the double poles located at Uj,Vj (j = 1 . . . N). Cala- 
brese and Cardy proposed the Ansatz 

(*,-(.„%)...*;(*,%», = {t[w))k '" = ™ {fn "' w] ' (2A1) 



5 



LI 

A 

k-l / \ k+1 

k \ 



u 

t 

V 

k-l J k+1 

Y 



Figure 1: Three branch points generated by the fusion of two oriented intervals 
(«i,fi) (1*2)^2) m the limit \u\ — u 2 \ — > 0. A general symmetry of the replica approach 
described in Section [3TT1 allows to deform the intervals as indicated on the right. 



with 

N , x 1 

; u j ^ u k , Vj ^ v k V j ^ k . (2.12) 

The condition Uj 7^ u k , Vj 7^ v k V j 7^ k, which is an important ingredient to derive 
(11.71) . has no justification in the process of sewing the n replicas, as there is no geometrical 
obstruction in letting two branch points of the same type, say u, coincide, like in Figure [TJ 

The equality illustrated in Figure [1] comes from the fact that, in the physical systems 
we are studying, a cut joining two branch points may be continuously deformed, as we shall 
discuss in detail in the next Section. Note that the cut associated with the fused branch point 
sews together the sheets in the order k — > k + 2. Clearly the order in which the replicas are 
sewn together is immaterial, provided that all the n replicas are cyclically connected. On the 
other hand, if n is an odd number, the permutation k — >• k + 2 covers once all the n replicas. 
Therefore, according to the general argument of [10], the complex scaling dimensions of the 
twist field arising from the fusion of two twist fields with equal charge are expected to 
be exactly the same, i.e. A n = A n — ^(1 — 4?). We can take advantage of this remarkable 
property to write down the operator product expansion (OPE) 

*+(0, 0) iD) = —^-S+^O, 0) + . . . (2.13) 

where n is an odd number. The dots indicate less singular terms. 
Similarly (12. 8p yields 

*+(O J O)*;(t O) t0) = r ^- + ... (2.14) 

Comparing (I2.13P and (I2.14p shows that two branch points of the same type develop, as 
the interdistance decreases, an UV singularity with an exponent which is exactly half the 
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exponent associated with a pair of branch points of opposite type, while fll .Tj) would predict 
a zero. A numerical check of this behaviour is illustrated in Figure [TIH 

In the limit configuration depicted in Figure HJ tr p n A becomes of course proportional to 
the n— th power of the three-point function 

trp n (u, v u v 2 ) oc ($; 2 (n, u), Vi), $+(w 2 , u 2 ))c > ( 2 - 15 ) 

therefore conformal invariance completely fixes its functional dependence: 



1 

\(vi - v 2 )(v 2 -u)(u- t>i)p 



trp n (u,v l ,v 2 ) oc — — — , (2.16) 



which is valid for any odd n. A numerical check of the scaling behaviour of this quantity for 
n = 3 can be found in Figure [91 

Going back to (12.111) and (12.121) . note that they are based on two distinct assumptions: 

i) (T(w))ii nN can be written as the Schwarzian derivative of a suitable function f n N] 

ii) this function is given by (12.121) . 

It is easy to prove that the latter assumption is not correct in the case of multi-intervals 
because the function f n ,N>i(w) does not fulfil the condition that the only singularities of 
its Schwarzian derivative {f n ,N,w} are the double poles at Uj and Vj, as conformal Ward 
identities require. A simple way to see it, is to put n = 1 in ( 12.121) . thereby eliminating all 
the branch points and the associated twist fields. If the only singularities of (12.1 ip were the 
double poles at Uj,Vj (j = 1 . . . N) one should have {f\N,w} = 0. This would imply, as 
emphasised at the beginning of this Section, f±N(w) to be a linear fractional transformation, 
which evidently is not the case when N > 1. Actually we find, in the N = 2 case, 

3 1 31 3/1 1 \ . _ 

{fi2,w} = --- - + , (2.17) 

2 (w — r\Y 2 [w — r 2 y r\—r 2 \w — r\ w — r 2 J 

where 7*1 and r 2 are the two roots of the equation = 0, namely 

v if 2 - u x u 2 ± a/(mi - v^h - u 2 )(u 2 - v 2 ){v 2 - m) 

r X2 — ■ . 

ui — v 1 + u 2 - v 2 

Since f n2 = (fi 2 Y^ n , it easy to verify that these double poles at are present, with the 
same coefficient, for any value of n. Therefore the regular zeros of f' n2 , where the mapping 
w f n2 is not invertible, behave like (fake) primary operators of complex scaling dimensions 
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S n = S n = —jq- Hence {f n 2,w} has unwanted singularities besides those associated to the 
branch points. 

On the contrary, in the case of a single interval it is straightforward to write down the 
inverse function w(z) = fnl( z ) which solves the ordinary differential equation (ODE) 

w' n oc {w-u) n ~\w-v) n+1 , (2.19) 

showing that the only points where this conformal mapping is not invertible are just the 
branch points u and v. 

Having proved that when N > 1 (12. 12ft is not a solution of (12.111) . the question arises 
as to whether there is any solution of such an equation. A crucial observation is that the 
n— sheeted Riemann surface TZ n ,N with N > 1 is topologically inequivalent to the complex 
plane C. A fortiori there is no way to conformally map TZ n .N>i onto C. As a consequence, 
one cannot longer use the main physical motivation outlined at the beginning of this Section, 
namely that the vanishing of (T)c leads to express {T(w))-ji n N>1 as a Schwarzian derivative. 
However this argument does not exclude a priori that (T(w))-ji n N>1 is a Schwarzian derivative 
for some other reason. Thus we can try to find possible solutions of (12. lip . 

From a mathematical point of view the problem can be reformulated as follows: find a 
function z(w) such that i) its Schwarzian derivative {z, w} is singular only at the branch 
points «i,M2, • • • ,U2N of R n ,N', H) these singularities are double poles; in) the coefficient^ of 
(w — Ui)~ 2 is 1 ~ 1 2 //rt . We shall see that this problem admits very few solutions which can be 
explicitly constructed. 

As a starting point of this analysis, notice that it is always possible to map a connected 
domain of C (with some boundary identified) onto R n ,N- For instance 7?.2,2, the double cover 
of the plane w branched over u 1; u 2 , m 3 , u 4 , is conformally equivalent to a torus, represented 
by a parallelogram T> of the z— plane with opposite sides identified. The ODE (I2.19p is now 
replaced by 

w' 2 oc (w — ui)(w — u 2 )(w — u 3 )(w — u 4 ) . (2.20) 
Its general solution may be brought into the form 

w(z) oc -^r + ui , (2.21) 

where p(z) is the elliptic p— function of Weierstrass associated with the mentioned paral- 
lelogram. Since the singularities of {z, w} coincide with the zeros of w' the first condition is 

2 The last condition implies that the scaling dimension of the associated primary field is A„. 
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fulfilled. An explicit, straightforward, calculation yields 



{z,w} 



4 

E 

i=l 



A in — ?; • ^— * 



(w — uA 2 Aw — Ui Ui — Uj 



(2.22) 



which is the sought after behaviour. It shows, once again, that in the TZ 2 , 2 case the solution 
of (jMID is not (12TT21 . but the inverse of the function (12~2TT) @. 

One may envisage a straightforward extension of (I2.19P and f)2.20p by looking for solutions 
of the equation 

w' n oc P(w) , (2.23) 

where P(w) is zero or singular only on the branch points Uj, vj (j = 1, . . . , N). It may be 
worth observing that evaluating {z, w} does not require the actual knowledge of the solution 
of (12.231) : using the identity {z, w} = —{w, z}/w' 2 we obtain 

(2ra-l)P 2 -2nPP 

y^- 2 , (2-24) 

where P = ^ and P = ^. 

aw dw 

Contrarily to naive expectations, there are very strong restrictions on the possible form 
of P(w) if we are to exclude branch points associated to primary fields with negative scaling 
dimensions. First, P(w) must be a polynomial of degree 2n at most. However, if the degree 
is less than In there is a branch point at infinity on the n— sheeted Riemann surface, thus 
we assume that P(w) is exactly of degree 2n in w, but this is not yet sufficient. It is worth 
noting that the absence of negative scaling dimensions coincide with the requirement that 
there are no movable singularities |15j . i.e. singularities of the solution w(z) which move in 
the z— plane as the initial values are varied. 

If we demand either the absence of movable singularities or the positivity of the scaling 



dimensions of the involved twist fie 
of allowed ODE of the kind (T2T231 



ds, it turns out that there are essentially only six types 
j. Two of them are (l2~T9l) and (l2~20ft . A third type is 
simply a limit case of (12.201) when u\ — U2- The other three types are 

u/ 3 oc (w -u^tw -u 2 ) 2 (w -u z ) 2 , (2.25) 
w' 4 oc (w — Mi) 3 (w — u 2 ) 3 (w — m 3 ) 2 , (2.26) 
w ,6 oc {w-u 1 ) 5 {w-u 2 ) 4 {w-u 3 ) 3 . (2.27) 



3 This inverse function may be expressed in terms of an elliptic integral of the first kind. 
4 See for instance [TS], p.312ff. 
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They describe multi-sheeted Riemann surfaces branched over three points as depicted in 
Figure [TJ In particular the first, once inserted in (12.241) and combined with the conformal 
Ward identities, leads to ( 12.161) with n = 3. The latter two correspond to n = 4 and n = 6 
and could therefore be used for hints to extend the OPE (12.131) to even n. 

Comparison of (12.221) with the conformal Ward identity associated to (T(w)ti 2 2 suggests 

tr p 2 (u 1 ,u 2 ,u 3 ,u 4 ) oc | — |e/12 ■ (2.28) 
j<k<4 \ u i Uk \ 

This expression has the expected symmetry properties under the permutations of the branch 
points; note also that rescaling all the coordinates by a common scale factor U{ — ► A Uj yields 
trp 2 (Au,j) = A~ c / 2 p 2 (uj), which is the expected scaling property of the square correlator of 
four twist fields (see also (14.51) ). However Eq. (12.281) cannot be trusted because it implies, as 
it would be simple to show, the unwanted identity (T(z))x> = 0, where T> is the parallelogram 
of the z— plane conformally equivalent to 7?.2,2- 

Summing up, we have shown that Eq. (12. lip has no other solutions besides the two-point 
and some three-point correlation functions. Thus the problem of finding the entanglement 
entropy of disjoint intervals remains open. 



3 tr p n A as a vacuum expectation value 

In this section we explicitly write tr p n A in whatever quantum system as the vacuum expec- 
tation value of a suitable observable O, defined on a larger system composed of n decoupled 
copies of the original system. This method is particularly well suited for numerical simula- 
tions because a single numerical experiment directly yields the value of trp^. 

The partition function Z = tie' 1311 of our d— dimensional quantum system at inverse 
temperature j3 can be computed in a standard way by doing the Euclidean functional integral 
in a d + 1-dimensional hyper-cubic lattice A = {x, r} (xj, reZ) over fields <f)(x) = 0(x, r) 
periodic under r — > r + (3. Therefore the system composed of n independent replicas of the 
original system is described by the n— th power of Z: 

/n 
\[V{<p k ]e-^ s ^ , (3.1) 
k=i 

where <pk is a field configuration associated with the k— th replica and S[4>] is the Euclidean 
action. Though the method has a much wider applicability, we assume, for the sake of 
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Figure 2: The branch points in a square lattice are located on the sites of the dual lattice. 
The links intersecting the cut A connect two consecutive planar lattices labelled by k and 
k + 1 (modn). 



simplicity, that the fields 0& are associated with the nodes of the lattice and that S is the 
sum of contributions of the nodes and the links of the lattice A 



S[ ( f )k ] = Y, v ^{x)) + Y,F{<p k {x)Ak{y)) 

(x y) 



(3.2) 



with (xy) ranging over pairs of adjacent nodes in A; V and F are arbitrary functions. 

To define the coupled action S^[<f)i, 02 • • • , 4>n], let us begin with a two-dimensional 
system defined in a square lattice. The subsystem A consists of one (or more) segments 
joining pairs of nodes in the dual lattice A. We now connect these pairs of nodes with a line 
A, not necessarily coinciding with the segments of A, and replace each link intersecting this 
line with a link connecting two consecutive sheets, k and k + 1, in a cyclical order (see Figure 
[2]). In this way we obtain a discretized version of a n— sheeted Riemann surface where the 
two points on the dual lattice are the branch points and the line A is the cut. The associated 
coupled action is 



Q (n) 



E 

k=l 



\ x )l 0fc+I(mod 



(3.3) 



The generalisation to higher dimensions is almost obvious: in a 2 + 1 dimensional system 
the two (or more) branch points are replaced by one (or more) closed paths 7 on the dual 
lattice and A is replaced by a surface E whose boundary is 7, and so on. 

The observable we are interested in is 

Its vacuum expectation value in the system of n independent copies of the original system is 



Z n {A) 



to Pa 



(3.5) 
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Z n (A) is the partition function of the system in the n— sheeted Riemann surface or of its 
multi-dimensional generalisations. The reduced density matrix pA is normalised as in [TO] . 

In numerical simulations this observable can be evaluated with great accuracy, because 
only the links intersecting the cut A (or its multi-dimensional generalisation) contribute to 
it. Actually, it is not necessary to simulate n independent systems: simulating a single 
system and taking, for each measurement of O, n statistically independent configurations is 
enough. Alternatively, as a consistency check of the algorithm, one may simulate the system 
on a n— sheeted Riemann surface (or its multi-dimensional generalisations) and verify that 
in such a case one should obtain 



3.1 A useful symmetry 



A system composed of n decoupled copies of the same d— dimensional quantum system 
has an interesting invariance which is something more than the standard replica symmetry 
considered in disordered systems. Let us assume that the system in its Euclidean description 
is defined on a stack {A& , k = 1, . . . n} of n copies of a d + 1 dimensional hyper-cubic lattice. 
Each node of the stack is characterised by d + 2 integral coordinates (xi, X2, ■ ■ ■ , Xd+ii k) with 
1 < k < n. Assume furthermore an imaginary time coordinate Xd+i running form — oo to oo 
(i.e. zero temperature) and define the following transformation 

s i - a ; i (t = l,2,...d+l) ; *-{£_ 1(modn) tZllta (3 ' 7) 

where a is an arbitrary real number. This transformation may be thought of as a simple 
relabelling of nodes of the stack: after such a transformation the system consists always of 
n disjoint lattices, with exactly the same geometric structure as before the transformation; 
the only difference is that the label k is no longer constant along a given copy. As a conse- 
quence, the partition function Z n of the composed system is obviously invariant under such 
a transformation. Note however that it is not a symmetry of the classical action S[(j)f.]: 
one has to perform the functional integration in order to implement this invariance; therefore 
in a sense it can be considered to be a quantum symmetry. The above considerations can 
be generalised in an obvious way to the coupled system described by the partition function 
Z n (A). 

To make the discussion concrete and explicit, we specialise to the case where the field 
theory in question is defined on an infinite cylinder of circumference L (see Figure [3]). The 
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Figure 3: The k labels of the nodes of the stack of n cylinders after the transformation (13. 7p . 



transformation just defined could be also generated by the n— sheeted Riemann surface 
described in Figure [2] by moving the two end points in opposite directions in such a way 
that the cut A does wrap around the cylinder. In the limit where the end points touch each 
other, they annihilate and disappear, and the n— sheeted Riemann surface becomes a stack 
of n disjoint cylinders. Thus, the symmetry generated by the transformation (13.71) can be 
viewed as the invariance of the system under the addition (or the removal) of closed cuts 
wrapped around the cylinder. We can enlarge the game by also considering cuts associated 
with homotopically trivial loops, which can be used for instance to show that the cut A 
connecting two branch points my undergo continuous deformations, or to prove identities of 
the kind 

trPA={(Mi,i>i) = t r PA={(ui,Vj)(uj,Vi)...} ■ ( 3 - 8 ) 

This point is illustrated in Figure HI In words, these identities show that the cuts or intervals 
associated with the branch points Uj and Vj are in no way distinguished lines on the surface: 
their introduction has a similar role as the choice of a reference frame on the surface. Note 
that the essential physics responsible for this symmetry is not some special effect of CFT, 
but rather a general property of the functional integration - a basic tenet of generic quantum 
field theories. 

If we combined together open and closed cuts associated with arbitrary permutations of 
replicas, a very rich structure would emerge. However it goes beyond the scope of this paper. 
In the following, we consider for convenience only cyclical or anti-cyclical permutations of 
the replicas. 

Two points u and v on a cylinder can be connected by a cut A in two topologically 
inequivalent ways (see Figure E]). These two ways correspond to the two complementary 
subsystems A and B of the whole quantum system under study. Now if we combine an 
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Figure 4: A graphical proof of the identity (13.81) . 




Figure 5: A graphical proof of the equality trp^ = trp^. 

open cut associated to the cyclic permutation k — > k + 1 with a closed cut associated with 
the anti-cyclic permutation k — > k — 1, like in Figure the mentioned symmetry evidently 
interchanges the role of the two subsystems A and B and yields at once the fundamental 
identity 

trp^ = trp£ ,Vn . (3.9) 

It is known that this relation, as well as its ensuing consequence (11.21) . is valid only if the 
whole system A U B is in a pure state. This equality is violated at finite temperature 
because a thermal state is in a mixed state, of course. As a consequence we can infer that 
the mentioned symmetry should break down at finite T. For, since the quantum system 
corresponds to a two-dimensional Euclidean theory with a compactified periodic imaginary 
time with period f3 = 1/T, the infinite cylinder reduces to a torus of size L x /3, and the 
transformation (13.71) is no longer a symmetry of the system, as it transforms n disjoint tori 
into a single torus of size L x n/3. 

In our numerical experiments we used the vanishing of the difference tr p\ — tr as a 
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sort of thermostat to keep the temperature of the simulated system low enough. 



4 Entanglement in Q state Potts Models 



One of the simplest and more studied models of statistical mechanics is the Q— state Potts 
model [T6l [T7] . which is the basic system symmetric under the permutation group of Q 
elements. It can be defined by associating with each node x of an arbitrary lattice A the 
field or spin variable <j> x — 1, 2, . . . , Q . Its partition function at the coupling J is taken to 
be 



where Sq[4>] = — J2( xy ) J ^<t>x<t> y with (xy) ranging over pairs of adjacent nodes on A. 

In a two-dimensional lattice this system has a typical order-disorder transition which is 
continuous in the range < Q < 4. Contrary to the naive expectation, the clusters made 
of adjacent sites with aligned spins do not play an important role at criticality. A different 
definition of cluster was proposed [18]. These clusters are defined as adjacent sites with the 
same spin connected by bonds with probability p = 1 — e~ J . Within such a definition, these 
clusters behave correctly at the critical point, in the sense that their radius and the density 
of the percolating cluster scale with the correct critical exponents. 

The partition function (14. ip can be rewritten in terms of these clusters using the Fortuin 
Kasteleyn representation [I~9] : 



where v = = e J — 1; the summation is over all spanning subgraphs of A, namely the 
subgraphs made with all the nodes of A; b(G) is the number of edges of G, called active 
links, and c(G) the number of connected components or Fortuin-Kasteleyn (FK) clusters. 
This formulation of the partition function, sometimes called di-chromatic polynomial, enables 
one to generalise Q from positive integers to real and complex values. In particular Q = 2 
corresponds to the Ising model and Q — 1 is the random percolation problem. 

The implementation of the general method described in Section [3] to evaluate tr p n A is very 
simple in this case, even if the action in the FK formulation has no longer the form (I3.2p . as 
c(G), the number of clusters, is a non-local quantity. Each configuration of the Q-state Potts 
model is uniquely characterised by the location of the active links. Consider a stack of n 
configurations of this kind, defined on a square lattice. Choose an interval A connecting two 
arbitrary nodes of the dual lattice and sew together the n sheets according the rules drawn 




(4.1) 




(4.2) 



GCA 
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in Figure [2] so as to form a covering of the n— sheeted Riemann surface with a connected 
square lattice. Note that the sewing operation does not modify the local structure of the 
lattice nor the location of the active links, while the number of clusters c(G) may change. 
Therefore, according to (13.51) . we obtain the explicit and simple result 



where Ck denotes the number of clusters of the fc-th replica before sewing and ca is the total 
number of clusters in the stack of n replicas sewn together along A; the vacuum expectation 
value is taken with respect the system composed of n decoupled copies of the Potts model 
under study. 

Although so far we have considered for the sake of simplicity a two-dimensional system, 
it is clear that the formula (14.31) still holds true in any space dimension, provided one defines 
appropriately the subsystem A. 

A first obvious consequence of ( 14. 31) is that in random percolation, i.e. Q = 1 Potts model, 
trp^ = 1 for any A and any n: as intuitively expected, there is no quantum entanglement 
in random percolation. 

Some further remarks are in order. The Q-state Potts model on a square lattice for 
continuous Q varying between and 4 undergoes a second order phase transition at v — \/Q. 
Its critical behaviour is described by a CFT with a conformal anomaly c related to Q by [20] 



Random percolation at the percolation threshold corresponds to c = 0, as the absence of 
entanglement combined with (I2.10p requires. 

For Q > 1 there are very efficient non local cluster Monte Carlo algorithms which can 
be applied for integer Q [2"Tl |2"2"] as well as for for continuous Q [23J. They are particularly 
well suited for accurate estimates of tr p\ through the formula (14.31) , as at the heart of these 
algorithms there is a reconstruction of the FK clusters of the configurations. 

We simulated in this way a two-dimensional critical Ising model, which corresponds to 
a CFT with c = |. We considered very elongated lattices of size L x L' with periodic 
boundary conditions on either side. Typically, the aspect ratios L'/L were between 4 and 
8. This choice amounts to a convenient compromise between (I2.10p . which is expected to 
be exactly true only on an infinitely long cylinder, and the numerical experiments, which 
are necessarily performed on a finite system. As a criterion to decide whether our cylinders 
were long enough, according with the discussion following (13.91) we compared the values of 



trp^ = {Q CA ~^ Ck ) 



(4.3) 




(4.4) 



1(3 




0.32 I 1 1 1 1 1 LJ 

4 6 8 10 12 14 16 

r 

Figure 6: Correlation function of two twist fields in a stack of n = 3 independent copies of a 
2D critical Ising model in a square lattice of size 32 x 256. The data are generated by 2 • 10 6 
Monte Carlo configurations. The solid line is a one-parameter fit to (I2.10p . 

trp^ and trp^ , where A and B are two complementary cuts along the circumference of 
the cylinder, as in Figure [51 When V / L = 8 we found tr p n A ~ tr p 1 ^ within the statistical 
errors. The behaviour of this quantity at n = 3 as a function of the distance of the branch 
points is depicted in Figure The solid line is a one-parameter fit to (I2.10p . where the only 
free parameter is the unknown proportionality constant. We conclude that our numerical 
experiments reproduce very well the expected conformal behaviour of the Tsallis entropy, 
even if the actual size of our lattices is not very large. 

When considering the entanglement entropy of critical systems corresponding to CFT 
with negative c, an intriguing aspect emerges: according to the general expression (11.31) . Sa 
becomes negative. One is led say that the mixed state obtained by tracing over the degrees 
of freedom of the complement of A is in some way more ordered than the vacuum \^>), 
which has vanishing entropy. While all that may sound counterintuitive, it has a simple, 
testable counterpart in the replica approach: according to (12.101) . the correlation function of 
the branch-point twist fields should grow with their distance. One can check it on critical 
Potts models with Q < 1, since they correspond to c < 0. The non-local cluster algorithms 
[2"Tj |2"2"| |2"3] are applicable only for Q > 1, however we can resort to a local algorithm [21] 
which can be implemented in an efficient way in the range < Q < 1 [25], [26] . 

We performed our numerical experiments at Q = 2 — y/3, whose critical behaviour is 
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Figure 7: Correlation function of two twist fields in a stack of n = 2 independent copies of 
a 2D critical Q— state Potts model, with Q = 2 — y/3, in a square lattice of size 25 x 100. 
The corresponding CFT has a negative conformal anomaly (c = —W). The data come from 
2 ■ 10 6 Monte Carlo configurations generated with the algorithm described in [26]. The solid 
line is a one-parameter fit to (12.101) . 

described by the (non-unitary) minimal model M 7;12 , which has c = — according to (I4.4p . 
Also in this case we found good agreement with the predicted behaviour of the two-point 
twist fields correlator, as Figure shows. 

One of the most general and fundamental properties of correlation functions of local 
operators at criticality is the scaling behaviour under a common rescaling of all the dimen- 
sionful quantities. We may use this principle to test one of the most important findings of 
the replica approach to entanglement [10] , namely the discovery that tr p n A is proportional to 
a correlation function of twist fields. Although the specific form of tr p n A in the case in which 
the whole system is enclosed in a square box of side L and the subsystem A consists of N 
disjoint intervals A = (ui, Vi), . . . , (un, Ujv) is not actually known, (11.61) implies the following 
homogeneity property (see Figure [8]) 

trp n (Xu 1 , Xvx, . . . , Xu N , \v N , XL) = X~ 4nAnN tr p n («i, v u . . . , u N , v N , L) , (4.5) 

where A is any positive rescaling factor and A n is the scaling dimension (11.51) of the 2N 
branch point twist fields. We checked this formula in critical Ising systems in the cases 
N = 1 and N = 2, finding perfect agreement for boxes large enough, as Figure [8] shows: 
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Figure 8: Comparison of the scaling behaviour of the twist correlator in the case of two and 
four branch points. The x axis is in a log scale and the numerical data lie on straight lines 
as expected. According to (14. 5p the slope is exactly 4nA n N, where 2N is the number of 
branch points, while the intercept is the only fitting parameter. In this instance the critical 
model is the Ising model and n — 2, therefore 4nA n = |. 



the slopes of the two straight lines match precisely with the predicted value 4nA n N. We 
checked in the same way the scaling properties of the three-point function (12.161) proposed 



in the present paper, which corresponds to put N = § in the above formula, finding again a 
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good agreement (see Figure [9]), even if such a case is more demanding from a computational 
point of view. 

We also performed a numerical experiment to settle whether in the limit |wi — W2I — > 
the observable tr p\ = u ui Vl \i U2 „ 2 )} vamsri es, as predicted by (11.71) . or diverges, as argued in 
Section [2j For this purpose we considered two collinear intervals, as sketched in Figure 
[TUl on a stack of 3 replicas. In order to reduce the noise, instead of considering as in the 
other numerical experiments three independent Ising systems, we simulated the model on a 
three-sheeted Riemann surface IZ3 with a cut in (u±,vi), thus the measured observable was 

w^wwk - m^^pmMh . (4 . 6) 

Figure fTUl shows that this quantity grows as \u\ — W2I decreases, even if within this set up 
it is not possible to accurately measure the associated exponent, which according to (12.131) 
should be exactly half that of a pair of branch points of opposite charge. 
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Figure 9: Scaling behaviour of the twist correlator in the case of three branch points of a 
critical Ising model defined on a stack of n = 3 replicas in the configuration depicted in 
FigHJ The x axis is in a log scale and the numerical data lie on a straight line as expected. 
The slope of the solid line is exactly 6nA n = ~ as predicted by (12.161) . The data of larger 
lattices are generated by ~ 10 7 Monte Carlo configurations. 

5 Conclusions 

In this paper we discussed both analytically as well as numerically some properties of the 
trace of the n— th power of the reduced density matrix p& in the special case in which A is 
a subsystem of a quantum system described by a conformal field theory in two dimensions. 
From the analytical point of view we pointed out that when A is composed with more 
than a single interval the explicit formulae proposed in [10] suffer from some inconsistencies. 
The origin of these can be traced to the fact that in the derivation of the general formulae 
certain singularities of the Schwarzian derivative have been overlooked. Of course a moot 
derivation does not imply a wrong final result, thus it is well possible that the formulae for 
the entanglement entropy obtained by analytic continuation to n — > 1 are still correct. 

From the computational point of view we developed a simple method which allows us 
to directly evaluate trp^ as a vacuum expectation value. We applied it to two different 
critical systems in two-dimensional Euclidean lattices, corresponding to two different values 
of conformal anomaly c, finding perfect agreement with the predicted formulae in the case 
in which A is composed by a single interval in a system of finite size [TO] . 
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Figure 10: In a critical Ising model defined on a stack of n = 3 replicas, the subsystem 
A is chosen to consist of two collinear segments (ui,Vi) (u2,v 2 ) as indicated in Figure and 
the quantity ( 14. 6 j) is plotted versus |tti — 1x2! ■ The size of the lattice is 128 x 128 and 
I iti — v\\ = 1 1*2 — v 2 \ = 32. The data are generated by 2. 10 6 Monte Carlo configurations. 



A distinguishing feature of the present method is that the simulations are performed in the 
unperturbed system: the n replicas of the system under study do not interact with each other 
nor are sewn together along some particular subsystem A. The information on entanglement 
is encoded in the observable O whose evaluation does not perturb the system. An obvious 
advantage of such an approach is that one can exploit the same set of configurations to 
obtain information on the entanglement entropy for a variety of subsystems. 

We would also like to remark that the usual method to measure c in lattice simulations 
is very different from that one may infer from the present paper. The standard method 
originates from the observation that the first subleading correction to the free energy of a 
CFT on a long cylinder is universal and proportional to c [H] , however the free energy cannot 
be directly measured in lattice simulations, thus usually one circumvents this difficulty by 
measuring and integrating the internal energy, with an inevitable loss of precision. The 
method described in this paper may be viewed as a novel and powerful tool to directly 
measure c: the rescaling property (14.51) of the Tsallis entropy tells us that a measurement of 
tr p n A in two lattices of different size suffice to fix the slope of the upper straight line drawn 
in Figure El and hence the value of c. 

To conclude, notice that the proposed numerical method can also easily be applied to 
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critical and non-critical systems in any space-time dimensions. For instance, adding one more 
dimension to the Ising system simulated in this paper and applying Kramers and Wannier 
duality [27], one may obtain information on the entanglement entropy of a confining gauge 
theory, an issue of growing interest in the last years [281 12H1 EQJ I3T] . 
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